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^USTR^CT 


have  ana]>zed  the  energy  conversion  processes  occurring  in  three  din'i'. 
sional  driven  recofinect  ion .  In  particular,  we  studied  the  energv  conversion 
processes  during  localized  reconnection  in  a  tail  like  rr^aqnet  ic  corf  i  gur  at  i  or, . 

'*e  find  that  3D  driven  reconnection  is  a  powerful  energv  converter  which 
transforms  magnetic  energy  into  plasma  bulk  flow  and  thermal  energv.  Three  di¬ 
mensional  driven  reconnection  is  an  even  more  powerful  energy  converter  than  ?1) 
reconnection  because  in  the  30  case  plasmas  were  drawn  into  the  reconnect  i  ori  re¬ 
gion  from  the  sides  as  well  as  from  the  top  and  bottom.  Field  aligned  currer,ts 
are  generated  by  30  driven  reconnection.  We  have  identified  the  physical 
mechanism  responsible  for  these  currents  which  flow  from  the  tail  toward  the 
ionosphere  on  the  dawn  side  of  the  reconnection  region  and  from  the  ionosphere 
towards  the  tail  on  the  dusk  side.  The  field  aligned  currents  form  as  the  neu¬ 
tral  sheet  current  is  diverted  through  the  slow  shocks  which  form  on  the  outer 
edge  of  the  reconnected  field  lines  (outer  edge  of  thf  plasma  sheet). 


I.  INTHODUCTION 


A  magnet ospher ic  substorm  is  a  good  example  of  an  explosive  energy  conver¬ 
sion  process  in  cosmical  plasmas.  Over  the  past  two  decades,  two  models  of  the 
substorm  have  evolved,  the  unloading  model  or  energy  storag'^  and  release  mode) 
and  the  driven  model.  In  the  unloading  model,  reconnection  on  the  dayside 
magnetopause  leads  to  enhanced  convection  and  the  storage  of  magnetic  energy  in 
the  magnetotail  magnetic  field  (McPherron  et  al . ,  1973;  Russell  and  McPherron. 
1973).  This  leads  to  thinning  of  the  near  earth  plasma  sheet  and  to  the 
formation  of  a  new  neutral  line  in  a  localized  region  in  the  near-earth  (10  ~  lb 
R^)  plasma  sheet  (McPherron  et  al . ,  1973;  Russell  and  McPherron,  1973;  Hones, 
1973).  This  reconnection  impulsively  releases  the  stored  energy  at  substorm 
onset  to  power  geomagnetic  activity.  The  driven  model,  too,  starts  with 
magnetic  merging  on  the  dayside  magnetopause.  However,  in  this  model  the 
magnetosphere  responds  directly  to  variations  in  the  external  solar  wind 
conditions  (Perreault  and  Akasofu,  1978;  Akasofu,  1979,  1950).  In  particular, 
the  energy  input  rate  from  the  solar  wind  equals  the  energy  output  rate  during 
the  substorm  in  the  driven  model. 

In  this  paper,  we  present  a  three-dimensional  model  of  tail  reconnection 
during  substorms.  Our  primary  purpose  is  to  s.udy  the  energy  conversion 
processes  occurring  in  the  tail.  When  we  wish  to  contruct  a  sound  physical 
model  to  explain  any  naturally  occurring  phenomenon,  it  is  essential  to 
understarid  the  fact  that  the  system  is  energetically  "open".  To  reveal  an 
energy  conversion  process,  the  first  step  we  take  is  to  extract  a  portion  (black¬ 
box)  out  of  the  causal  chain  in  which  the  process  we  are  concerned  with  is 
taking  place.  The  black  box  must  be  extracted  so  that  there  is  little  feedback 
of  energy  from  the  output  to  Input  during  the  course  of  the  energy  conversion. 
Otherwise,  the  model,  no  matter  how  mathematically  sound,  is  piiysically 


unjcctpt  db  1  f- .  If  tbf  bldPk  box.  is  so  sr^ioil  tfidt  d  X  d' t  i  d  1  p.irt  of  t  Ot 

no' I'ut  f-nf  rg\  relurns  to  l  tio  input  on  t  bo  t  i  mr  scdlf  of  iUrrest,  then  ttif  rr  -o' 

ohxiodsK  iririuclfs  3',  intern,)!  i  ricc)ns  i  st  f  nrv  .  Tht-  n>  >1  r  .st  t.)e  expdooe:!  sc 
t  hdt  t  t)e  output  energy  returning  to  the  input  Odo  be  inrluctep.  Or.  I  tie  ottier 
hand,  if  the  black  box  is  so  large  that  several  separablr  eriergy  conversior. 
prc'cesses  are  included,  then  the  model  is  accept  able.  The  global  simulatiori  of 

tlie  rr.agnet  osphere  is  an  example  of  suet'  a  system  arid  in  principle  is  the 

preferred  v^ay  to  study  the  substorm  process.  Honsever,  because  of  the 
limitations  of  the  available  computers,  the  model  is  often  so  coarse -ora i nec 
that  the  phvsical  process  underlying  the  energy  conversion  can  be  overlooked. 

In  order  to  comprehend  the  physical  process,  therefore,  it  is  desirable  to 
divide  this  black  box  into  self-consistent,  smaller  boxes. 

In  a  substorm,  a  substantial  part  of  the  energy  output  is  lost  in  ttie 

iemosphere.  The  lost  energy  goes  into  the  heating  of  neutral  atoms  and 

moleculf-s,  so  that  it  is  unlikely  that  ttie  lost  energy  goes  back  to  the 

magnetosphere  during  the  course  of  a  suhstorm.  Ori  ttie  ott.er  hand,  simil.ir 
s.il)Slorns  are  observed  repeatedly.  From  the  repeatability  of  substorms 

t ‘,f  existence  of  definite  energy  loss  in  the  ionosphere,  we  ca^  cor.cluce 
I 'idt  ir.ere  must  exist  a  finite  energy  input  to  the  system.  Accordingly,  a  model 
s,:l''.t  cT.m  must  be  constructed  to  include  a  proper  iriput  and  output.  From  ca.jsal- 
itv,  ttie  input  must  change  ttie  conditions  wittiin  the  black  oox  i  ri  sucti  a  w.iy 
t  uat  ,)  sutistorin  is  t  r  i  (|<i(  red  . 

iv<  li.ive  desi'’|n<(l  our  irx.d)  I  to  investigate  ttie  effect'  ot  co'. vert  ion  ori  the 
1  "  I  !  i  fi :  i  on  of  tail  reco.niK'Ct  i  .oil .  In  ttiis  ciriven  r  e-ronriec  t  i  fm  mode- 1  , 
re  I'li  .ff'iion  results  from  flow  from  ttie  lolie-s  normal  tu  the  neuiral  stie  e  I  (bato, 
l'.-‘7','i.  Ill  order  to  h.fve  sufficie-nl  resolotion  iri  l  tu  m  ale  1  ,  we  consifleT  ori  1  v  ,i 
S's.i!!  volume  around  ttie  recnnne'Ct  ion  re-gion.  The-  input  from  ttie  solar  w  i  nr!  is 
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determined  by  the  boundary  conditions  at  the  top  and  bottom  of  the  simulatin' 
box.  The  term  driven  (»hen  applied  to  reconnection  does  not  have  the  same 
meaning  as  <»hen  it  is  used  to  describe  phenomenological  mc'dcls  of  substorms.  ^ 
»»e  Hill  see  in  this  study,  driven  reconnection  is  an  active  energy  converter  in 
which  magnetic  energy  stored  in  the  tail  lobes  is  converted  into  plasma  enerqv . 
Thus  the  driven  reconnection  model  is  similar  to  the  unloading  model  of 
substorms . 


IJ  .  ue  UFfi-\L  HEMARkS  OU  FMLKG'i  COtULHSJOU 

()(-riving  the  enrrqs  f^qualions  from  the  ordiriary  set  of  resistive  mdqoeto- 

tivdrodvnamir  (MHD)  equations,  we  obtain 


pv  •  dS  +  I  3  X  B  •  _^dV  -  j  Vp  •  ^dV 


jl  =  -  J  •  dS  +  /  7p  •  vdV  +  / 

,  E  X  B 

-  f  = - =  •  dS  -  /  3  X  B  •  vdV  - 

where  K  =  /  -1  pv^d\'  (kinetic  energy) 


r,32d; 

(2) 

/  n3^dV 

(3) 

(4) 

(thermal  energy) 


(3) 


W  =  f  —  dV  (magnetic  energy)  (6) 

and  dS  is  the  surface  element  vector  and  dV  is  the  volume  element.  In  (1)-(6), 
p  is  the  mass  density,  v  is  the  velocity,  3  is  the  current  density,  B  is  the 
magnet  ic  field,  p  is  the  pressure,  y  is  the  polytropic  index,  r,  is  the  resistiv¬ 
ity,  and  E  is  the  electric  field. 

These  equations  give  us  all  the  necessary  information  about  the  macroscopic 
energy  conversion  in  a  black  box.  The  left-hand  sides  of  (1)-(3)  represent  ttie 
changes  in  the  total  kinetic,  thermal  and  magnetic  energies  in  the  black  box. 

The  first  terms  on  the  right  represent  the  output  and  input  powers  of  the  black 
box.  As  described  in  the  introduction,  t  tie  input  power  is  essential  in  the 
present  model  since  it  is  responsible  for  tfie  change  in  ttie  initial  equilibrium 
configuration  which  triggers  an  energy  conversion  process  in  the  black  box.  The 
output  power  provides  a  key  to  the  energy  conversion  process  because  we  are 
interested  in  how  much  energy  is  ejected  from  the  black  box  as  a  consequence  of 
ttie  conversion. 
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The  last  two  terms  on  the  right-hand  sides  of  (1)-(3)  represent  the  mutual 
energy  exchanges  among  the  kinetic,  thermal  and  magnetic  energies  inside  the 
black  box.  They  can  give  us  information  about  the  efficiency  of  the  energy 
conversion  process.  They  can  be  used  to  answer  questions  such  as  how  quickly 
the  magnetic  energy  is  converted  into  kinetic  energy.  As  a  consequence  of  this 
internal  conversion,  the  Internal  energy  is  ejected  out  of  the  black  box.  The 
ejected  (output)  energy  can  in  turn  be  an  input  energy  for  the  succeeding  event 
in  the  external  causal  chain.  For  instance  in  the  magnet ospheric  substorm,  this 
succeeding  event  can  be  the  excitation  of  auroras. 

A  diagram  of  energy  conversion  processes  is  shown  in  Figure  1. 


Ill  .  LKTLHtaLL'i  DKULtj  KECUUIIL XT  ION 

In  the-  recormpction  process  the  energv  source  to  be  released  is  the  tr.aqnel - 
ic  energy.  In  the  magnet ot a i I ,  the  magnetic  energy  consists  of  two  components. 
One  is  the  dipole  field  and  the  other  is  the  field  produced  by  the  diamagrietic 
current  flowing  in  the  plasma  sheet.  The  dipole  field  is  a  potential  field  and 
cannot  be  released  at  all.  The  only  free  energy  is  the  magnetic  energy  stored 
in  ttie  diamagnetic  currents  which  deform  the  dipole  to  form  the  tail.  Thus  an 
external  energy  agent  such  as  the  solar  wind  which  deforms  the  magnetosphere 
must  play  a  leading  role  in  triggering  reconnection  in  ttie  plasma  sheet. 

It  would  be  natural,  therefore,  to  presume  that  the  solar  wind  intruding 
into  the  magnetosphere  through  dayside  reconnection  could  drive  reconnection  in 
the  magnelotail.  The  intruding  solar  wind  effects  the  tail  in  two  ways.  First, 
the  magnetic  free  energy  is  stored  in  the  plasma  sheet.  This  is  manifested  as  a 
plasma  sheet  thinning.  Second,  the  invading  plasma  triggers  reconnection  by 
compr.’ssing  the  plasma  sheet.  To  approximate  this  process,  we  developed  the 
driven  reconnection  model  [Hayastii  and  Sato,  197S;  Sato  et  al.,  1978;  Sato  and 
Hayashi,  1979;  Sato,  1979;  Sato  et  al . ,  1983]. 

Indeed,  this  idea  that  reconnection  would  occur  somewhere  in  the  magneto¬ 


tail  has  been  demonstrated  by  global  simulations  of  the  solar  wind-magnetosphere 
interaction  [Leboeuf  et  al.,  1981;  Brecht  et  al.,  1983;  Wu,  1983]. 


U  .  [jUlLiUCAL  MOULL 

The  numerical  mode)  has  been  described  previously  [Sato  el  al.,  1S'^3:  S.ii  o 
and  Havashi,  1962].  Here  we  will  summarize  il  briefly  for  convenience. 
started  with  a  very  simple  equilibrium  configuration.  The  magnetic  field  and 
plasma  configurations  we  used  are 

B(x,y,2)  =  (Bplanhiz/L ) ,0,0) 

—  Pq 

T(x,y,z)  =  Tosech^iz/L) 

where  B  is  the  magnetic  field,  p  is  the  mass  density,  T  is  the  temperature.  H^,. 
Dq,  Tq  are  constants  and  L  is  the  half-width  of  the  neutral  sheet  in  the  z 
(north-south)  direction. 

He  assumed  that  reconnection  is  triggered  by  a  non-uniform  compression  of 
the  plasma  sheet.  In  order  to  compress  the  plasma  sheet,  we  injected  a 
magnetized  plasma  through  the  north  and  south  boundary  planes  (z  =  t  L^)  of  the 
simulation  box.  Initially,  the  incoming  mass  flow  pattern  was  such  that 

pv2(z=+L2)  =  +-j^(1  +  cos  ^)(1  +  cos  j^) 

4  L,^  Ly 

where  is  the  z-component  of  the  velocity,  L^,  Ly,  art  the  dimensions  of 
the  simulation  box  and  is  a  constant,  pv^^  and  pv^  were  initially  set  to 
zero.  Later  they  were  adjusted  so  that  the  mass  flow  vector  was  always  perpen¬ 
dicular  to  li  al  the  boundary.  The  other  bound.iries  (x  -  t  and  v  -  L.  )  were 
open  boundaries  through  wliich  plasmas  can  freely  enter  or  exit.  He  reduced  ttm 
physical  domain  of  the  simulation  box  by  assuming  symmetry  (or  ant i - sy mmet ry 
when  appropriate)  about  the  equator  (z=0),  tl»e  noon-midnight  meridian  (yrO)  and 
tt>e  neutral  line  (x=0).  He  have  adopted  the  resistivity  (n)  model  [Sato.  1979] 

n  =  a(Vy-  Vf.)^  for  Vp  >  V^, 


m 


where  Vp  is  the  diamagnetic  drift  velocity  and  a  and  are  constants.  Vq  r 
3/ne  where  3  is  the  neutral  sheet  current,  n  is  the  plasma  densitv  and  e  is  ttn 
electron  charge. 

The  simulation  system  was  a  rectangular  box  with  dimension  L^=3L,  L^=5L  and 

=  2L .  This  was  implemented  on  a  41  x  40  x  point  grid. 

In  the  actual  calculat ions,  all  variables  were  normalized  to  the  following 

Bo 

parameters:  L,  Length;  =  - - ri/2»  velocity;  Bq,  magnetic  field;  Pq,  mass 

tUoPo’ 

density;  Bq/uqL,  current;  pressure;  B^^/uq,  energy;  YqLV;^,  resistivity; 

and  Bq\ electric  field.  In  the  computer  run  presented  here  we  set  Aq=.2,  a- 
0.02  and  V(.=3  while  Bq,  Tq,  Pq  and  Vq  were  all  normalized  to  1. 
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V  .  IIIJMIKICAI  WtSULTS 

In  the  tl)  model,  only  t  tie  compress  i  onal  maqnelohy  drod>  nam  ic  modes,  ttie  fast 
arid  slo*s  modes.  Can  be  taken  into  account,  but  the  sfiear  Alfvt'n  mod<  canriot  . 
Consequently,  no  field  aligned  currents  can  be  modeled  because  it  is  the  Mfvfn 
*»ave  that  can  carry  the  field  aligned  current.  In  order  to  model  field  aligned 
Currents,  a  3D  model  is  necessary. 

That  field  aligned  currents  can  be  generated  by  3D  externally  driven 
reconnection  has  been  demonstrated  b>  Sato  et  al  .  (19S3).  Starting  «»ith  the  MHl) 
momentum  equation,  the  continuity  of  current,  mass  continuity  and  the  MHD 
approximation,  the  field  aliqried  currents  in  the  magnetosphere  can  be  calculated 
b>  solving 


cs'b  ’  B 


VB  -  3 


1  n 


(b*t7)v  + 


(7) 


w  i  t  li 


Pi  =  V  X  ^ 

i  i  .2in 


and 


=  cB  X  VP/B' 


X  (£,^1 


wh(  re  3.  and  3j^  are  the  parallel  and  perpendicular  current  densities,  B  is  the 
maqrietic  field,  p  is  the  mass  density,  iJ,  is  t  tie  parallel  component  of  P  the 
vorticity,  is  the  perpendicular  component  of  '.i,  ^  is  tiie  inertia  current, 
ft  is  ttie  density,  P  is  the  pressure  and  v  is  ttie  velocity  (Haseqaiva  and  Sato, 
1979;  Sato  and  lijima,  1979;  Vasyliunas,  19S4).  The  partial  derivative  on  t  ti< 
left  tiand  side  of  equation  (7)  is  alonq  H.  Sato  (  1982)  has  argued  ttiat  ttie 
inertia  current  term  is  negligible  and  that  the  agent  of  the  firld  aligned 
current  must  be  the  3^  •  VB  term  in  equation  (7).  The  diagram  stiowing  his  pre¬ 
diction  is  reproduced  in  Figure  2. 
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He  drgupd  lhal  ds  reconnfct  i  on  occurs  iti  a  locdiized  region  of  th(  tail, 
two  slow  shocks  will  form  into  which  the  cross  tail  current  will  be  diverted. 

In  addition,  B  will  increase  at  the  slow  shocks  due  to  tne  pile  up  of  fielc 
liries  which  are  driven  towards  the  neutral  sheet.  Thus,  •  VB  >  0  in  tne  dawri 
sector  and  •  VB  <  0  in  the  dusk  sector.  The  resulting  field  aligned  currents 
are  directed  from  the  tail  towards  the  ionosphere  on  the  dawn  side  and  outward 
from  the  ionosphere  towards  the  tail  at  dusk. 

In  this  paper,  we  will  use  the  simulation  results  to  verify  ttiis  predic- 
tinr..  Contours  of  the  generated  field  aligned  current  intensity  in  two  y-z 
cross  sections  (x=2.23L  and  xz3.0L)  are  shown  in  Figures  3  and  4.  For  reference 
reconnection  occurs  at  x=0.  Two  pairs  of  field  aligned  currents  are  generated; 
one  is  in  the  northern  hemisphere  (z>0)  and  the  other  is  in  the  southern  hemis¬ 
phere  (z<0).  On  the  dawn  side  (y<0),  the  field  aligned  current  flows  from  the 
magnetosphere  towards  the  ionosphere,  while  on  the  dusk  side  it  flows  from  the 
ionosphere  towards  the  magnetosphere  as  predicted  by  ttie  model.  Two  additional 
features  should  be  noted  in  these  figures.  First,  the  field  aligned  currerd 
d^nsitv  is  as  large  as  ttie  initial  neutral  sheet  current  derisity  in  a  small  re¬ 
gion.  This  indicates  ttiat  locally  the  neutral  sheet  current  is  almost  complete¬ 
ly  interrupted.  Secondly,  the  peak  positions  of  the  field  aligned  currents 
expand  in  t  tie  z-direction  as  the  observation  point  moves  along  the  x-axis  away 
from  the  x  line.  This  suggests  ttiat  the  field  aligned  current  is  generated  at 
the  slow  shock  layer  [Sato  et  al.,  lyh3]. 

In  order  to  prove  that  the  field  aligned  current  originates  from  ttir-  J,  • 

''.tl  term,  we  tiave  plotted  the  neutral  stieet  current  3j  superimposed  on  ttu-  con¬ 
tours  of  ttie  magnitude  of  Itie  magnetic  field  B  in  Figure  S.  If  we  examine  ttie 
current  distribution  witfi  care  we  find  that  tlie  current  derisity  arrows  in  t  tu- 
eguatorial  plane  (z=0)  vanisti  at  the  edges  of  the  reconnection  region. 


F  urt  fif-rmorc ,  thf  diverqenrf  of  tfo-  rurrf'ots  tircurs  arross  lf>f  B  contours.  If. 

[  i()urr  6  v*tif-rc  have  overlapped  Figures  4  and  b,  we  cari  clearlv  see  that  ttie 
co'lours  of  the  field  aligned  current  density  coincide  with  the  di\erqence  of 
the  neutral  sheet  current.  The  field  aligned  currents  occur  where  3j^  •  7B  #  0 
as  predicted  by  Sato's  model.  Note  that  the  field  aligned  current  density  is 
verv  small  (J,,  <  .02)  in  the  equatorial  plane  in  the  region  where  the  cross  tail 
current  reverses  direction.  The  field  aligned  currents  are  limited  to  the 
region  of  the  slow  shocks. 

lAe  also  have  investigated  the  contributions  to  from  the  terms  in  (7) 
wich  are  proportional  to  the  vorticity,  Q.  The  parallel  vorticity  term  gives  a 
field  aligned  current  which  is  opposite  to  that  observed.  Therefore,  term  (3) 
is  not  important  for  generating  these  currents.  Both  terms  involving  the 
perpendicular  vorticity  are  appoximately  zero  and  are  much  smaller  than  the 
first  term.  Thus  we  can  conclude  that  the  field  aligned  currents  are  generated 
b\  the  3,  •  VB  term. 

-  X 

In  equilibrium  the  pressure  contours  are  parallel  to  the  magnetic  field 
intensity  contours.  This  means  that  the  diamagnetic  current,  which  is 
perpendicular  to  Vp,  must  flow  along  an  equi-  B  contour  in  an  equilibrium. 

Thus,  the  field  aligned  current  generation  is  due  to  the  non-equilibrium 
magriel ic  field  configuration  caused  by  externally  driven  reconnection. 


As  noted  above,  externally  driven  reconnection  produces  slow  shocks  in 
wtiich  high-speed  plasma  flows  are  found  [Sato,  1979].  In  thf-  21)  simulation,  ttu 
flows  reached  the  Alfvfn  velocity.  The  3D  simulation  confirms  a  very  similar 
bftiavior  of  plasma  flow  at  least  in  the  noon-midnight  meridian  plane  [see  f igurf 
2  of  Sato  ef  al.,  1983].  In  3D,  the  plasma  flow  can  become  super-magnet oson i c . 
In  Figure  7,  we  have  plotted  the  flow  speed  distribution  (V,^)  along  the  x-a\is. 


Thf  curve  Vp  represenls  lt>f-  dislribuHori  of  ttir  fast  mugru  l  obofi  i  c  speeri .  It- 
plasma  flo*v  is  almost  linearly  accelerated  along  the  x.-avis  from  the  reconnec¬ 
tion  region  (x=0)  up  to  about  x=1.  Since  the  velocitv  is  normalized  by  the  inj. 
tial  Alfven  speed  in  the  upstream  region,  ne  see  that  the  speed  well  exceeds  ttip 
Alfvfn  speed  beyond  x  ^  1.  The  flow  speed  exceeds  the  local  fast  magnetosonic 
speed  (V'p)  beyond  a  certain  point  (x  >  x^)  and  the  local  magnetosonic  Mach 
number  (V,^/Vp)  reaches  approximately  2.  The  local  Mach  number  of  the  flow 
increases  monotonica 1 1 y  with  distance  beyond  x^.  However,  the  absolute  flow 
speed  reaches  a  maximum  at  the  super-magnet oson ic  point  (x::x^)  and  ther, 
decreases.  This  implies  that  some  deceleration  mechanism  is  operating  beyond 
that  point.  The  neutral  sheet  current  density  (3^)  distribution  is  plotted 
along  the  x-axls  in  Figure  6.  For  x  >  2,  the  current  reverses  direction  (see 
also  Figure  5  of  Sato  et  al.,  (1963)).  Since  the  reconnected  field  has  a 
positive  (northward)  z-component  in  the  positive  x  region,  the  _2  force 
points  towards  the  x-point  in  the  reversed  current  region  and  the  plasma  is 
decelerated  in  this  region.  This  caused  the  velocity  decrease  in  Figure  7. 

Sato  et  al.  (19S3),  examined  the  curren,ts  in  the  equatorial  plane  and  found  a 
large  current  vortex  between  the  reconnection  region  and  the  region  of  reconnec¬ 
ted  field  lines.  The  current  pattern  in  Figure  7  is  along  the  symmetry  axis  of 
this  vortex. 

Let  us  now  examine  the  corresponding  distributions  for  the  pressure  (P)  and 
the  reconnected  field  .  There  are  two  very  important  features  in  Figure  9. 
Mrst,  the  reconnected  field  increases  linearly  with  distance  for  x  £  xq.  Ttiis 
is  a  clear  indication  of  the  diffusion  region  predicted  by  the  steady  reconnec¬ 
tion  theory  [Petschek,  1964;  Vasyliunas,  1976].  Note  that  the  plasma  is  almost 


linearly  accelerated  as  a  function  of  x  in  this  region  (see  Figure  7).  Ttie 
theory  further  predicts  that  beyond  the  diffusion  region  both  the  reconnected 
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fif  1(1  r)(i(J  tti(  prfssure  slm'jlrl  be  curislard  in  a  rfyiofi  ra  1  1  f-ri  t'a  rt'()i(it.. 

Siirh  a  pliJtfau  regio-'i  is  observed  in  tbr  simulation  for  xq  ■  \  ^  x^.  Most 

S'irpr  i  s  i  ng  1  \  ,  however ,  is  .i  large,  sharp  df-erease  in  t  hf  prt  ss-jre  pri^file  at  x  •- 
Xj.  This  is  the  point  where  t  tie  local  magnet  oson  ic  Mach  number  exceeds  one 
(see,  Figure  7).  A  close  examination  of  the  profile  reveals  that  when  ttie 
sharp  nrop  in  the  pressure  occurs,  also  starts  decreasing.  Recall  that  the 
flow  velocity  decreases  for  x  >  x^  (Figure  7).  The  steep  decrease  in  pressure 
suggests  that  this  may  be  a  fast  shock  propagating  in  the  positive  x  direction. 
However,  there  are  difficulties  with  this  interpretation.  The  density,  which  we 
would  expect  to  decrease  if  this  is  a  shock,  increases  -  *^9).  In 

addition,  although  the  magnetic  field  decreases  as  we  would  expect,  the  decrease 
is  very  small  -  1-03). 

F.nergv  Conversion 

In  this  section,  we  examine  the  energetics  of  3D  driven  reconnection.  Let 
us  first  examine  the  relationships  between  the  input  and  output  powers  which  are 
aiver,  by  t tie  first  terms  on  the  right-hand  side  of  (1)-(3).  Since  there  are 
three  independent  boundary  planes,  the  output  boundary  iy-z),  the  side  boundary 
(x-z)  and  the  input  boundary  (x-y),  we  have  9  input  and  output  power  values  to 
consider.  These  9  values  are  plotted  as  a  function  of  time  in  Figure  10.  The 
three  curves  with  the  first  three  numbers,  1,  2  and  3,  correspond  to  the  kinetic 
energy.  The  three  curves  with  the  second  three  numbers,  4,  S  and  6,  correspond 
to  the  ttiermal  energy.  The  remaining  three  curves,  7,  S  and  9,  correspond  to 
ttie  magnetic  energy.  Ttit  smallest  number  in  each  group,  1,  4  or  7,  represfriis 
t  tic  output  power  through  t  tif  output  boundary.  The  middle  number  in  each  gruup, 

,  “>  or  S,  represents  ttu  output  power  ttirouqli  the  side  boundary.  The  last 
number,  3,  6  or  9,  represents  the  input  power  ttirougti  Itif  input  boundary. 


A  positive  value  for  output  po»%er  indicates  that  eriergy  flows  out  throuqti  thi 
bouridary,  while  a  negative  value  indicates  that  energy  flows  in.  A  positive 
value  for  the  input  power  represents  energy  flow  into  ttie  sif^ulation  reaion  and 
a  negative  value  represents  outflow.  F rrxn  the  input  curves,  3,  6  and  9,  we 
notice  that  most  of  the  input  power  is  carried  by  ttie  magnetic  field.  Ttiis  may 
well  represent  the  actual  magnet ot ail  where  the  energy  is  accumulated  predomi¬ 
nantly  in  terms  of  the  magrietic  energy. 

The  input  powers  in  3,  6  and  9,  stimulate  reconnection  and  energy  conver¬ 
sion  in  the  black  box.  The  resulting  output  kinetic,  thermal  and  magnet ic 
powers  are  shown  by  1,  4  and  7.  Since  the  induced  super-magnet osonic  flow  is 
limited  to  the  region  within  the  two  slow  shock  layers  (i.e.  the  model  plasma 
sheet),  where  the  thermal  energy  is  predominant,  the  dominant  output  power  is 
the  thermal  power.  The  first  rise  in  the  thermal  power  (4)  (0  <  t  <  8t;^)  is 
caused  by  the  convection  system  set  up  in  the  black  box  from  the  input  bouridary. 
The  second  enhancement  in  the  output  power  starting  at  T  ~  is  induced  l)y 

reconnection.  Coincident  with  this  second  rise  in  the  thermal  power,  ttie  kine¬ 

tic  power  is  greatly  increased.  In  contrast,  the  change  in  ttie  output  magnetic 
power  is  not  so  remarkable.  Since  the  output  power  caused  by  reconnection  is 
limited  to  the  slow  shock  regions,  the  calculation  of  the  output  power  for  ttiis 
discussion  also  should  be  limited  to  ttie  slow  shock  region.  However,  this  does 
not  change  our  conclusions  since  over  80%  of  the  output  power  is  wittiin  the  slow 
shocks.  Clearly,  the*  total  output  power  caused  by  reconnection  exceeds  ttu 
total  input  power.  Thus,  the  energy  conversion  process  occurring  in  the  tilack 

box  is  an  active  process,  not  a  simple  pumpitui  process  induced  by  the  input 

flow.  Now  let's  look  at  the  energy  outputs  ttirougti  the  side  t)oundary,  ?,  8  arid 
b.  There  is  very  little  kinetic  energy  input  or  output  (?).  However,  srimf 
power  is  pumped  out  of  the  black  box  througti  the  side  boundary  early  in  ttie 


oimi.ildl  iod  tiefor*-  rtTu'id!  rt  io'i  hjs  t)«'qufi .  Durincj  tin  r  <  ('udit-r  I  ii)i(  slaiic,  tin 
pi'vvcr  ttiroijijti  t  fip  side  douridarv  beromes  nf=-qative.  This  mf-dds  that  thir  s  i  dr 
tio'jddar)  has  brromr  an  rriergs  sourer  durinq  rrronnrot  i  (>'■  .  T  ri  rlarifv  this  we 
presprit  the  flow  pdltrrn  in  the  y-z  plane  at  x=()  (Tr14.6)  in  f  iqure  1.  Reroi'- 
nert ion  is  occurring  in  the  center  of  the  figure.  It  is  evident  that  plasma  is 
be  i  no  sucked  into  the  reconnection  region.  Ttiis  gives  additional  strong  evi¬ 
dence  ttial  externally  driven  reconnection  is  a  powerful  eriergv  civnverler.  I  r, 
t  tie  previous  2D  simulations,  no  suet,  ptipnomenon  was  alloweci.  This  is  a  ttiree- 
dimensional  effect.  It  is  this  inflow  through  the  side  boiridary  which  makes 
driven  reconnection  in  t  tiree-d  i  mens  i  ons  a  more  powerful  er.prgy  converter  than 
tt.at  in  two-dimensions.  The  fast  stiock-like  structure  observed  in  Figure  9  may 
be  attributed  to  this  effect  of  three-dimensionality. 

In  order  to  study  the  energy  conversion,  we  need  to  examine  the  work  done 
within  the  simulation  box.  There  are  tt.rec  sources  of  work,  the  work  done  by 
ttie  3  X  B  force,  the  work  done  by  tlie  pressure  force  and  the  ohmic  heating. 

Ttipse  works  are  plotted  in  Figure  1?.  In  the  early  pumping  or  convection  si  aip' 
(0  <  t  <  Si/y),  both  the  3  x  B  force  and  the  pressure  force  behave  almost  identi¬ 
cally,  which  means  that  no  net  acceleration  of  plasma  occurs.  The  oscillations 
are  caused  by  the  bouncing  of  the  magnetosonic  wave  between  the  input  bound¬ 
aries.  Tture  is  a  marked  differ^-nce  in  the  curves  during  the  reconnection  st  agr 
(t  >  TOt/^).  The  net  difference  goes  into  the  acceleration  of  plasma  at  ttie 
expense  of  the  magnetic  energv  .  Ttie  small  positive  contribution  of  the  pressure 
work  implies  that  some  adiabatic  compression  is  occurring.  In  addition,  ohmic 
t. eating  contributes  to  the  heating  of  ttie  plasma  during  ttie  reconnection  period. 

Tt.e  temporal  behavior  of  ttie  total  kinetic,  therm.il  and  magnetic  energies 
are  stiown  in  Figure  13.  iarly  in  ttie  simulation  as  ttie  plasma  flows  from  ttie 
input  boundary  toward  the  neutral  sheet,  ttie  magnetic  energy  and  the  thermal 


Al  th(  samp  lime,  the  magnet  ir  enerqv  deviates  from  t^e  linear  qrowtti.  Tt,( 
differences  of  the  magnet  ic  and  thermal  energies  fro^'  the  initial  linear  gr(«:n 
line  give  the  total  amounts  of  energv  ejected  from  t'le  black  bov.  Ncitf  ttat  tft 
kinetic  energy  starts  to  increase  drastically  in  the  reconnection  stage. 


VI.  DlSCUbSlOIJ  AIID  CONCLUSlOUb 


Hri^fd  o'l  tiir,  ions,  S.ito  (  IV??)  rtrqufd  tb.d  dri\f'i  i c  r<  ri'''- 

ri>-'r'tio^i  could  b('  d  powerful  criC'rqy  coriverter  wfiich  t rd'tsf orned  mdj'ielic  enerav 
irito  pldsmd  bulk  flow  and  therndl  energy.  Such  an  explosive  energv  releasf 
meebdoism  is  necessary  to  understand  maqnt  t  ospher  i  c  substorir.s  or  solar  flares. 

In  this  study,  we  tiave  analyzf-d  t  be  energy  conversion  processes  occurring 
in  t  tiree-d  i  meris  i  ona  1  driven  reconnect  i  on  .  In  particular,  wr  investigated  ttie 
case  in  wtiicti  reconriec  t  i  ori  occurred  in  a  localized  region  in  a  t  a  i  1  -  1  i  ke  m.tgof  t  - 
in  (\)nfi(juta!  iofi.  Such  localized  reconrie-ct  i  ori  is  Itiougfil  t  rj  oci'ur  in  ttif-  m,3(jiie- 
totail  durina  substorms.  We  found  ttiat  30  driven  recorinect  i  ori  is  an  eve-ri  more 
powerful  energy  converter  than  the  20  case.  This  occurred  because  in  ttie  30 
e:ase  plasmas  were  drawn  into  ttie  reconnection  region  from  the  sides  as  well  as 
trie  top  arid  bottom.  In  the  30  case,  the  flow  velocity  reactied  the  magnetosonic 
velocity.  It  is  interesting  to  note  that  the  largest  velocities  occurred  at  the 
m.up'iet  f'soci  ic  point  witti  lower  velocities  fartlier  from  the  reconnection  region. 
rt<  found  that  the  neutral  stieet  currerit  reversed  signs  in  t  tif  regiori  of  redun  ti 
velocity  so  that  the  d  x  B  force  opposed  the  flow  and  dicelerated  ttu'  plasma. 

As  noted  in  the  i  n  t  rodur  t  i  on ,  the  term  driven  wtien  applit'd  to  reconnect  inn 
does  not  torn  the  same  meaning  as  when  it  is  used  to  descritic  pfa  riomeno  1  og  i  ca  1 
rriod'ls  of  sitistorms.  Driven  reconnec-l  i  ori  is  an  active  energv  cor. verier  iri  wliicti 
magrietic  energy  is  converted  into  plasma  tticrmal  energv  ,ind  flow.  In  P  igure  U), 
the  rate  of  at  wfiich  energy  is  released  is  murti  faster  t  tian  t  tir  rate  at  wtiicfi 
energy  is  input  into  the  model.  Thos  drivrn  reconnection  is  similar  to  ttu 
vuloadiiK;  model  of  svibsturms  in  wtiicti  magnetic  energy  stnrerl  in  ttie  tail  lobes 
is  converted  into  plasma  enerejy  . 

AnottifT  important  feature  of  ttie  driven  reconnee  t  i  nri  mode]  is  t  tie  fnrm.it  inn 
of  slow  stiocks  near  the  outer  edge  of  ttie  rrcnnnect  ed  field  li'a's. 


1  ta-  neut  r  ,1  I 


shtet.  curreni  is  diverted  Ihrouyti  ttif  slow  shocks.  This  leads  to  the  general  iofi 
of  field  aligned  currents  by  the  •  VB  term  in  equation  (7).  Sato  et  a). 

(  1963)  and  Walker  and  Sato  (  1953)  have  compared  tt/e  resulting  current  svst-^m 
with  that  inferred  from  observations  of  substorms  in  the  tail.  They  found  that 
the  current  system  in  Figures  3  and  4  with  field  aligned  currents  directed  from 
the  tail  toward  the  earth  on  the  dawn  side  and  from  the  earth  toward  the  tail  at 
dusk  is  consistent  with  the  current  wedge  model  most  frequently  used  to  inter¬ 
pret  substorm  observations  obtained  near  midnight. 

We  should  note  that  Birn  and  Hones  (1981)  have  done  a  similar  calculation 
and  find  field  aligned  currents  flowing  in  the  opposite  direction  to  those  in 
Figures  3  and  4.  In  this  model  reconnection  occurs  following  the  initiation  of 
anomalous  resistivity  in  in  the  tall.  Their  nwdel  does  not  include  the  initial 
convection  system  characteristic  of  driven  reconnection.  These  two  models  are 
not  necessarily  contradictory  and  it  is  possible  that  processes  analogous  to 
both  operate  in  the  tail  during  substorms.  For  instance,  driven  reconnection 
may  initiate  the  substorm  and  then  evolve  irtto  a  process  analogous  to  the  Birn 
and  Hones  model.  The  two  models  are  compared  in  greater  detail  in  Walker  and 
Sato  (1984). 
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r iQure  Capt ions 

Figure  1.  Schemut  ic  drawinq  of  tde  energy  conversion  processes  vs*"' i  c^'  ('.r.  'c.r 
in  driven  reconnection. 

Figure  2.  Schematic  drawing  of  the  plasma  sheet  dLirinq  r«conrirct  i  on  siuiw;-  ;  tt 
generation  of  field  aligned  currents. 

Figure  3.  Field  aligned  current  density  contours  in  the  y-z  (dawn-dusk)  plane 
at  X-2.2SL  at  1^14. 6T,y.  Ttie  values  oti  Itie  coiitoijrs  are  normalized  t  i.i  Ihf 
initial  neutral  sheet  currcfit  . 

Figure  4.  Same  as  Figure  3  for  X-3.0L. 

Figure  5.  Current  density  vectors  in  the  y-z  plaiip  at  xr3.0L  at  T-14.6i^. 

Superimposed  on  the  current  density  vectors  are  contours  of  magnetic  field 
magnitude. 

Figure  6.  The  field  aligned  current  density  contours  in  Figure  4  (stiaded 

region)  have  been  superimposed  on  the  current  density  vectors  and  maonetic 
field  contours  of  Figure  S.  Tlie  numbers  refer  to  values  of  the  field 
aligned  current  density  normalized  to  the  initial  neutral  sheet  curr(  d  . 

Figure  7.  Velocity  along  the  x-axis  at  T=14.6t/\.  Reconnection  is  occurrii>g  at 
x=0.  The  magrietosonic  velocity  is  given  by  . 

Figure  8.  The  cross  tail  current  density  along  tta  x-avis  at  Tr14.6';^. 

1  icjure  9.  The  pressure  and  rc’conruct  ed  field  aloni)  ttie  x-axis  at  1  =  14. 

figure  10.  Power  tlirough  tt\('  tlir<e  txmnd.ir  i  es  of  ttie  simulation  box.  1.  2  .iml 
3  give  the  kitud  ic  power  ttirouijl*  ttie  output  (x  tl  ),  side  (y-SL)  and  input 
(/=2L)  boundaries.  4,  *>  and  u  (|ive  the  thermal  power  ttiruu(]ti  ttie  tii.  ind.i- 
ries  while  7,  K  and  9  give  tlie  magnet  ie  power  t  liruugli  ttm  tioundarie'. 

Figure  11.  Flow  vectors  in  the  y-z  pi. me  at  x  =  0  at  1=14. 6t4. 


F  iqurr  1?.  V^ork  done'  v^ithjr,  the*  simulation  box.  ii'J  x  H)  •  \dV  give's  the*  evork 


done*  by  the  3  x  B  force*,  / VP  •  vdV  gives  the  **ork  done*  by  the  pressure 
force  and  ' '•i3'  dv  gives  the*  work  eione  by  the  3oule  healing. 

g2 

Figure  13.  Total  kinetic  (|p  v^dV),  thermal  magnetic  ^ 

energies  versus  time. 
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Abstract 


One  of  the  challenging  magnetospheric  plasma  physics  problem  is  the 
heating  of  heavy  ions  observed  in  the  equatorial  magnetosphere  in 
conjunction  with  ion  cyclotron  waves  (ICW)  generated  by  anisotropic  hot 
protons  (i.e..  GEOS-1  and  2  and  ATS  6  results).  The  complexity  of  the 
mechanisms  involved  in  this  phenomenon  (linear,  quas i - 1 i near .  nonlinear)  is 
such  that  a  numerical  simulation  is  the  only  way  to  disentangle  the 
different  plasma  characteristics.  We  present  here  the  results  of  such  a 
simulation.  The  plasma  parameters  which  have  been  chosen  are  those  which 
prevail  in  the  dayside  magnetosphere  at  geostationary  altitudes.  The 
Plasma  consists  of  cold  and  isotropic  and  He^  ions  with  a  small  number 
of  hot  anisotropic  protons  which  provide  the  free  energy  necessary  to 
generate  the  waves.  The  code  is  one-dimensional  in  length  and  three- 
dimensional  in  velocity.  It  is  electromagnetic  and  hybrid,  i.e..  the 
electrons  are  treated  as  a  massless  fluid.  The  results  obtained  during  the 
linear  phase  (t  <  150  )  are  in  agreement  with  those  expected  from  the 
linear  and  quasi -linear  theories  as  far  as  the  growth  rate  of  the  wave  (and 
the  frequency  of  the  most  amplified  wave)  and  the  variation  of  the  hot 
proton  anisotropy  are  concerned.  The  saturation  (B ~  0.^05  is 
explained  by  trapping  of  the  helium  particles.  But  the  most  interesting 
results  concern  the  heating  of  cold  species.  He"^  ions  are  heated  mainly  in 
the  perpendicular  direction  (kT  /2  150  eV.  kT,,/2  ==  80  eV)  and  they  are 

1  I 

heated  more  than  cold  ions  (kT  /2  ==  15  eV.  kT„/2  20  eV).  The  heating 
of  He*  ions  IS  a  two-step  process  :  first.  He*  ions  are  set  into 
oscillations  by  the  growing  wave  (in  both  v^  and  v^)  until  some  of  the  ions 
reach  a  parallel  velocity  of  the  order  of  the  resonant  velocity  at  which 
time  strong  heating  occurs.  Phase-space  plots  for  the  different  particle 
species  at  different  times  illustrate  the  time  evolution  of  this  heating 
mechanism. 
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It  has  recently  been  recognized  that  heavy  ions  are  important  to  the 
dynamics  and  stability  of  different  regions  in  the  earth's  magnetosphere. 
Shelley  et  al .  [1976]  and  Ghielmetti  et  al.  [1978]  have  reported 
observations  of  intense  fluxes  of  0^  streaming  up  high-lati tude  auroral 
field  lines  at  altitudes  of  ~  !  R^..  Data  from  the  S3-3  satellite  show  that 
ion  composition  can  vary  from  lOX  to  more  than  90%  oxygen  [Mizera  et  al . . 
1980]. 

The  inferences  of  0^  streams  deep  in  the  magnetotail  boundary  layer  by 
Frank  et  al .  [1977]  and  Hardy  et  al.  [1977]  suggest  that  there  might  be 
measurable  fluxes  of  energetic  0^  ions  of  terrestrial  origin  deep  in  the 
plasmasheet.  Results  from  the  energetic  ions  mass  spectrometer  on  ISEE-1 
show  that  the  plasmasheet  has  a  significant  and  variable  ionospheric 
component  (He"^  and  0^)  which  represents  from  more  that  10%  to  50%  of  the 
total  number  density  [Petersen  et  al..  19811. 

At  the  geostationary  orbit.  0*  and  He*  ions  of  medium  energy  i  -  10 
keV)  have  been  found  to  be  important  constituents  of  magnetospher i c  plasma, 
especially  during  magnetic  storms  [Geiss  et  at..  1978:  Young.  1979; 
Baisiger  et  al..  19801.  But  the  most  intriguing  result  is  the 
observations,  made  onboard  GEOS-1  and  -2  and  ATS-6  spacecrafts  that  He*  or 
0^  ions  are  heated  up  to  suprathermai  energies  (  -  100  eV)  at  times  when 
large  amplitude  ions  cyclotron  waves  (ICV's)  are  detected  [Young  et  al . . 
1981:  Mauk  et  al..  1981;  Roux  et  al..  1982;  Fraser.  19821. 

The  facts  that  heavy  ions  occur  in  many  different  regions  of  the 
magnetosphere,  that  they  are  of  terrestrial  origin  and  that  they  can  be  a 
major  constituent  of  the  plasma  motivated  the  community  to  reexamine 
theoretical  studies  of  wave  particle  interactions  including  the  effects  of 
heavy  ions.  When  studying  wave  particle  interactions  ir,  the  presence  of 
heavy  ions  there  are  two  problems  that  should  be  addressed.  The  first  is 
the  effect  of  the  heavy  ions  on  the  dynamics  and  stability  of  a  plasma. 


The  second  and  more  pertinent  problem  is  to  understand  the  mechanism 
responsible  for  energizing  such  heavy  ions. 

Several  studies  have  addressed  these  questions  for  different  regions 

of  the  magnetosphere.  On  auroral  field  lines,  the  heating  of  0*  has  been 

studied  by  using  electrostatic  waves  CLysak  et  al..  1980;  Papadopoulos  et 

al . .  1980;  Ashour-Abdalla  et  ai..  1981;  Ashour-Abdal la  and  Okuda.  19841. 

Lysak  et  al.  [19801  used  strong  turbulence  theory,  which  considered  the  ion 

orbits  in  a  set  of  fixed  amplitude  waves  with  a  definite  phase  relation. 

Papadopoulos  et  al.  [19801  argued  that  the  heavy  ions  form  a  minority 

constituent  of  many  magnetospheric  multi -ion  plasmas.  As  such,  the  heavy 

ions  do  not  affect  the  collective  mode  structure  of  the  plasma  due  to  their 

small  abundance  ratios.  They  thus  examined  the  acceleration  of  large  M/Q 

ions  in  the  presence  of  a  coherent  large  amplitude  electrostatic  hydrogen 

cyclotron  wave.  It  was  found  that  when  the  wave  amplitude  exceeds  a 

certain  critical  value  the  particle  orbits  become  stochastic  and  the 

particles  can  be  accelerated  by  the  wave.  The  most  important  result  was 

5/3 

that  the  maximum  energy  achieved  by  an  ion  scales  as  .  where  is 

the  mass  of  hydrogen  and  M  is  the  mass  of  the  ion  under  consideration. 
They  concluded  that  heavy  ions  are  thus  preferentially  accelerated. 
Ashour-Abdalla  and  Okuda  [19841  studied  the  acceleration  of  ions  on  auroral 
field  lines  associated  with  the  current-driven  electrostatic  ion  cyclotron 
waves  in  a  plasma  consisting  of  hydrogen  and  oxygen  ions  by  simulation 
techniques.  To  model  the  ionospheric  source  of  electrons  they  developed  a 
simulation  model  in  which  drifting  electrons  were  allowed  to  enter  at  the 
end  of  the  system.  They  found  that  the  oxygen  transverse  heating  generally 
exceeds  that  of  hydrogen  ions  at  the  low  altitude  end  of  the  field  line. 

There  have  also  been  studies  using  electromagnetic  waves  to  explain 
the  acceleration  of  heavy  ions.  Theoretically  it  is  well  known  that  the 
presence  of  heavy  ions  greatly  modifies  the  propagation  and  amplification 
characteristics  of  magnetospheric  electromagnetic  waves.  A  detailed 


analysis  of  the  wave  spectra  and  energetic  proton  distribution  nieasured  on 
GEOS-1  have  shown  that  the  observed  wave  spectra  can  be  interpreted  in  the 
framework  of  the  linear  theory  of  ion  cyclotron  instability  excited  by  an 
anisotropic  proton  distribution  in  a  plasma  containing  a  small  fraction  ~5 
-15%  of  He*  ions  [Roux  et  al..  19821.  Gomberoff  and  Cuperman  [19821  have 
also  performed  a  theoretical  analysis  in  order  to  evaluate  the  change  in 
the  temporal  growth  rate  of  such  instabilities  for  various  plasma 
parameters.  Gomberoff  and  “^eira  [19831  undertook  the  task  of  evaluating 
the  corresponding  change  in  the  spatial  growth  rate.  However,  such  studies 
have  generally  been  done  assuming  the  oackground  plasma  to  be  completely 
cold.  T(H*)  =  T(He+)  =  0.  whereas  it  ;s  known  from  'GEOS  and  DE  experiments 
that  both  electrons  and  ions  can  reach  a  temperature  equal  to  or  higher 
than  10  eV  [Decreau  et  al..  1982;  Chappell.  19831.  Lnaer  such  conditions, 
cyclotron  absorption  by  the  heavy  ion  species  will  compete  with  the 
amplification  induced  by  hot  protons,  and  growth  rates  will  be  strongly 
modified,  especially  i.n  the  vicinity  of  the  heavy  ion  gyrofrequency . 
Recently  Gendrin  et  al.  [19841  carried  out  a  parameter  search  of  ion 
cyclotron  waves  including  the  fiP’ie  temperature  effects  associated  with 
the  cold  species. 

The  heating  of  heavy  ions  by  electromagnetic  ion  cyclotron  waves  has 
been  addressed  oy  several  authors.  Gendrin  and  Roux  11980]  have  studied 
the  possibility  that  the  observed  heating  of  He*  ions  up  to  suprathermal 
energies  could  be  explained  by  quasi-linear  diffusion.  indeed  quasi- 
linear  diffusion  of  resonant  He’^  lOns  having  low  initial  parallel  and 
perpendicular  velocities  can  lead  to  a  considerable  increase  of  the 
perpendicular  He*  ion  distribution.  However,  to  heat  the  bulk  of  He* 
distribution,  sufficient  wave  energy  must  be  available  at  frequencies  where 
a  resonant  interaction  with  the  ccld  heavy  ion  species  can  take  place, 
i.e..  at  frequencies  which  are  near  the  heavy  ion  gyro: requency .  Gendrin 
et  al.  [19841  found  that  the  maximum  growth  rates  ''ccur  at  frequencies  far 


from  tne  heavy  ion  gyr:.f''equency;  ror.sequenily  thcie  waives  cannot  resonate 
with  the  bulk  of  heavy  ion  distribution.  In  search  of  another  heating 
mechanism.  Mauk  [1982.  19831  studied  the  nonresonant  interaction  between  an 
electromagnetic  ion  cyclotron  wave  and  thermal  He^  ions.  He  showed  that 
associated  with  the  increase  in  the  He*  ion  perpendicular  velocity,  phase 
bunching  was  observed.  Unfortunately.  Mauk's  calculations  were  done  by 
neglecting  the  ion  parallel  velocity,  which  led  to  incorrect  quantitative 
results:  an  apparently  infinite  increase  of  v  as  time  elapses  and  a  strong 
effect  even  when  the  wave  frequency  is  equal  to  the  helium  cyclotron 
frequency.  Berchem  et  al . .  [19831  and  Berchem  and  Gendrin  [19841  extended 
Mauk's  work  to  include  the  parallel  velocity  of  the  ions.  They  found  that 
phase  bunching  occurs,  and  the  particle  attains  large  perpendicular 
velocities.  They  carried  out  a  parameter  search  as  a  function  of  initial 
parameters  and  were  able  to  derive  approximate  analytical  expressions  for 
the  maximum  temperature  of  the  heavy  ions  attained.  These  calculations 
were  test  particle  calculations  in  which  the  wave  amplitude  was  kept 
constant,  and  were  therefore  not  self-consistent.  In  search  of  self- 
consistency.  Tanaka  and  Goodrich  [19841  used  simulation  techniques  to  study 
the  heating  of  heavy  ions  at  different  regions  of  the  magnetosphere  with 
particular  emphasis  on  the  heating  of  Alpha  particles  at  the  bow  shock.  As 
such  they  considered  a  plasma  consisting  of  hot  highly  anisotropic  protons 
in  the  presence  of  a  minority  species.  Since  these  studies  were  concerned 
with  gaining  an  understanding  of  He**  at  the  shock  region,  they  did  not 
include  a  cold  hydrogen  background  representative  of  the  ionospheric 
source. 

For  sake  of  completeness,  it  is  worth  mentioning  that  a  similar 
mechanism  (resonant  interaction  of  one  ion  species  with  waves  generated  by 
a  majority  ion)  has  been  invoked  to  explain  the  excess  of  He**  ion  bulk 
velocity  in  the  solar  wind.  However,  it  is  only  until  recently  that  the 
fundamental  modification  brought  in  to  the  dispersion  relation  by  the 


presence  of  the  heavy  ions  has  been  correctly  taken 


into  account,  which 


means  that  the  results  previously  obtained  are  not  relevant  to  the  problem 
we  are  studying  here  (see  Isenberg.  [1984a. bl  and  references  therein). 

In  this  paper,  our  aim  is  to  explain  the  heating  of  He*  observed  by 
the  (ZEOS  satellite  at  L  ~  7.  Our  starting  points  is  the  observation  set 
from  the  GEOS  satellite.  Having  developed  the  linear  theory  analysis  in  a 
recent  paper  [Gendrin  et  al . .  1984].  we  undertake  a  simulation  study  of 
electromagnetic  ion  cyclotron  waves  in  a  plasma  consisting  of  hot 
anisotropic  protons,  a  dominant  thermal  hydrogen  plasma,  and  a  minority 
species  of  helium  ions.  Following  this  introduction,  in  section  II  the 
computational  model  is  described.  simulation  results  of  our  three  ion 
component  plasma  are  presented  in  section  III.  We  conclude  in  section  iv 
by  summarizing  and  discussing  the  pertinent  heating  mechanisms. 


1 1  Simulation  Model 

We  have  developed  a  one  dimensional  electromagnetic  hyorid  simulation 
code  [Sgro  and  Nielson.  1976:  Byers  et  al..  I978J  where  electrons  are 
treated  as  a  massless  fluid  and  ions  are  treated  as  particles  moving  in  the 
four  dimensional  phase  space  (x.  v^,  v^..  v^).  Since  we  are  studying  low 
frequency  phenomena,  we  neglect  the  transverse  displacement  current  in 
Maxwell's  equations  (Darwin  approximation)  as 


J  =  rot  B  .  ( 2 ; 

We  also  neglect  electron  inertia  effects  for  low  frequency  waves,  and  we 
have  from  the  electron  momentum  equation 

-en^iE  +  "Jg  ^  3'  "  Po  =  '3) 

where  -e.  n  .  u  and  d  are  electron  charge,  number  density,  flow  velocity 
e  e  e 

and  pressure,  respectively.  Furthermore,  we  assume  the  quasi -neutrality 
condi 1 1  on 


-ene  *  q^n^  =0  '4. 

*nere  q  and  n  are  the  charge  and  number  density  of  an  "s"  ion  species, 

5  s 

respectively.  The  electron  pressure  is  integrated  in  time  by  using  the 
electron  energy  equation 


It  '  u„  ■  9rad)  =  -  I  Pe  “  “e 


Ihe  electron  flow  j  ;s  ootained  from  the  current  equation  as 

e 


u  =  (2q  njj 

w  S  5  S 


J )  /  n  q 
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inere  L'^  is  the  drift  velocity  of  a  "s"  ion  species.  The  ion  density  n^ 


and  drift  velocity  u^  are  calculated  from  the  motion  tf  ion  particles  whose 
PO.S  1 1 1  ons  and  velocities  are  determined  by  integrating  the  equations  cf 


mot  1  on 


dv 

dx 

dt 

_ 

"s 

1  E  +  V  1 

<  3  ). 

dt 

The  magnetic  field  B  and  the  electron  pressure  are  integrated  in  time 
using  (I)  and  (5).  while  the  electric  field  is  determined  by  a  predictor- 
corrector  method  using  (2),  13).  (4)  and  (6). 

For  simplicity,  spatial  variations  are  restricted  to  one  dimension, 
i.e..  the  X  direction,  and  periodic  boundary  conditions  are  assumed.  We 
focus  our  attention  on  waves  propagating  parallel  to  the  external  magnetic 
field.  Therefore,  both  the  wave  vector  k  and  the  static  magnetic  field  3^ 
are  in  the  x  direction.  We  consider  three  species  of  ions.  i.e..  cold  H*. 
cold  .He'*^  and  hot  with  bi -Maxwel  1 1 an  velocity  distribution.  Initially 
both  cold  components  are  isotropic  and  hot  protons  have  a  temperature 
anisotropy  of  T  'T  =2.  These  particles  are  distributed  uniformly  over  the 
simulation  space  consisting  of  64  grid  points.  No  wave  field  is  assumed 
initially  except  for  noises  introduced  by  particle  thermal  fluctuations. 

The  parameters  of  the  simulation  have  been  chosen  to  represent  typical 
conditions  prevailing  during  the  day  hours  at  the  geostationary  orbit  ^Tnere 
most  of  events  involving  lCW''s  and  energized  He*  ions  have  been  observed: 

-  magneti'c  field  intensity  :  140  nT. 

♦  .O 

-  coid  H  ion  density  !0  cm  '. 

-  coid  He*  ion  density  :  : 

-  thermal  energy  of  both  cold  species  :  l.TeV. 

-  anisotropy  of  both  cold  species  :  T  /T,,  =  !. 

-  hot  proton  density  :  I.!  cm”’ . 

-  thermal  parallel  energy  of  hot  protons  :  17  keV. 

-  anisotropy  of  hot  protons  :  T  /T,  =  2. 


Note  that  with  these  parameter  the  ratio  of  tne  :on  piasma  frequency 

to  the  proton  gyrof requency  is  equal  to  'OU.  and  tnat  the  Alfven 

■) 

energy  mV  “/2  -  5  KeV.  so  that  the  parallel  8  of  ’he  '-c,  a.sma  is  ^  0.3. 


In  the  simulation  we  ao  not  use  the  physical  parameters  as  expressed 

in  CGS  or  MKS  units,  but  instead  use  conventional  normalizations.  It  is 

simpler  to  use  a  system  in  which  the  Alfven  velocity  the  proton 

cyclotron  frequency  Q^.  the  proton  charge- to-mass  ratio  (q/ml^^  and  the 

magnetic  permeability  are  equal  1.  Other  parameters  are  normalized  to 

the  quantities  stated  above.  Since  the  static  magnetic  field  is  given  by 

2 

=  Qu/lq/min.  =  1  and  the  static  magnetic  energy  density  /2u  = 

0  n  rt  0  0  0 

0.5.  The  light  speed  c  is  equal  to  which  the 

•> 

electric  permittivity  -  i/c"  is  deduced.  The  thermal  velocity  of  the 
cold  protons  and  helium  are  0.018  and  0.009  respectively.  The  thermal 
velocity  of  the  hot  protons  is  1.8  and  2.55  for  the  parallel  and 
perpendicular  direction  respectively. 

The  grid  spacing  is  Ax  =  1  and  the  time  step  At  is  set  equal  to  0.05. 
The  number  of  grid  points  is  64  and  the  system  length  is  =  64.  The 
number  of  superparticles  in  the  system  is  4096  for  each  of  the  cold 
species.  As  for  hot  protons,  which  are  a  free  energy  source  in  the  system, 
a  relatively  large  number  of  superparticles  (32768)  is  necessary  in  order 
to  make  the  initial  thermal  fluctuation  noise  low. 
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The  GEOS  data  suggest  that  the  anisotropic  proton  distribution  excites 
ion  cyclotron  waves,  which  are  responsible  for  heating  the  thermal  helium 
ions.  In  order  to  test  this  hypothesis,  we  have  run  the  simulation  code  up 
to  Qj^t  =  1200.  We  will  first  consider  the  initial  phase  of  the  process. 

In  Figure  I.  we  see  the  time  history  of  the  magnetic  and  electric 

field  energy  densities  for  an  early  stage  of  the  simulation  run.  Both  the 

magnetic  field  and  electric  field  energy  densities  increase  throughout  the 

time  presented  in  this  plot.  The  perpendicular  component  of  the  electric 

field  IS  much  larger  than  the  parallel  one.  indicating  that  charge 

neutrality  is  conserved  by  the  predictor-corrector  code  and  that  the  wave 

IS  purely  electromagnetic.  It  should  be  noted  that  the  magnetic  field 

amplitude  attains  a  very  large  amplitude:  B  /E  =  0.055  at  CLt  =  200. 

wave  0  M 

Figure  2a.  shows  the  hot  proton  energy  density  in  a  direction  perpendicular 
and  parallel  to  the  ambient  magnetic  field.  Initially  at  time  a^t  =  100. 

the  temperature  anisotropy  decreases  only  slightly  to  (T  /T,).  =  1.9. 
Later  on.  the  temperature  anisotropy  decreases  more  dramatically,  as  the 
waves  attain  large  amplitudes,  so  that  at  Q^t  =  200.  (T  /T  =  1.6.  thus, 
the  hot  anisotropic  proton  distribution  is  responsible  for  the  growth  of 
these  waves.  The  lower  panels  (2b)  and  (2c)  show  the  energy  density  of 
thermal  helium  and  protons  respectively.  Both  thermal  ion  species  are 
accelerated  mainly  m  a  direction  perpendicular  to  the  ambient  magnetic 
field,  with  helium  ions  being  preferentially  accelerated.  We  will  discuss 
this  acceleration  phenomenon  later  and  we  will  show  that  it  is  in  fact  a 
two-step  process. 

To  identify  the  wave  mode  during  the  early  stages  of  the  simulation 
run.  we  plot  the  dispersion  characteristics  of  the  wave.  This  is  done  by 
Fourier  transforming  individual  k  modes  in  time.  Figure  3  is  a  plot  of 
frequency  w  versus  wavenumber  k  for  the  B^.  component  of  the  magnetic  field 
during  the  time  interval  Ct^t  =  1-204.  The  amplitude  of  the  wave  is  denoted 


by  the  height  of  the  peak.  It  is  well  known  from  linear  theory  analysis 
le.g.  Gendrin  and  Roux,  I98Q;  Gendrin.  1981)  that  the  addition  of  thermal 
helium  breaks  the  dispersion  relation  into  two  branches,  one  below  the 
helium  gyrofrequency  known  as  the  low  frequency  branch  (LF)  and  one  above 
the  helium  gyrofrequency  referred  to  as  the  high  frequency  branch  (HFK  In 
fact  the  LF  branch  and  HF  branch  are  clearly  depicted  here.  The  fastest 
growing  frequency  occurs  in  the  LF  branch  with  w/D,^  =  0.16.  =0.29 

and  B  /B  =  0.023.  The  dashed  line  is  the  real  part  of  the  dispersion 

w  o 

relation  calculated  from  the  linear  theory  for  the  parameters  at  Qj^t  = 
100.  By  comparing  the  dashed  line  and  the  plot  of  the  simulation  results, 
the  agreement  between  the  real  part  of  the  linear  dispersion  relation  and 
the  simulation  results  are  obvious. 

To  understand  the  time  evolution  of  the  wave  properties  during  the 
linear  stage  we  plot  the  time  history  of  four  dominant  modes  in 

Figure  (4a). 

The  value  of  k.  the  wavenumber  is  related  to  the  mode  number  m.  In 
fact  k  =  (27i/L^)ra  where  L^  is  the  system  length.  Thus  for  m  =  I.. 

=0.1  wnereas  for  m=2.  =0.2.  since  L^  =  64.  =  1  and  =  I.O. 

In  Figure  (4b)  we  plot  the  linear  growth  rate  versus  k  for  both  the  LF 
and  HF  branch.  From  the  upper  panel  we  see  that  the  mode  4. 
attains  a  large  amplitude  at  an  early  time.  The  growth  rates  of  these 
modes  are  found  by  calculating  the  slope  of  the  straight  dashed  line  shown 
in  Figure  4.  It  is  found  that  both  of  the  growth  rates  of  mode  3  and  4  are 
r 'Q.^  =  0.02.  Linear  theory  (Figure  4b)  predicts  the  peak  growth  rate  to  be 
r/Q^  =  0.022  and  to  occur  at  -  0-3  m  agreement  with  the 

simulation. 

Having  ensured  that  we  can  reproduce  and  understand  the  linear  theory 

from  our  simulation  study  we  now  examine  the  non-linear  physics  of  the 

interaction. 

In  Figure  5  we  plot  the  time  history  of  the  magnetic  and  electric 


hroughout  the  run.  As  noted  previously,  the  magnetic 


energy  densities  t 

field  increases  until  about  Cl,t  =  250.  reaching  a  peak  value  B  /B  =0.06 

n  wo 

and  then  remains  at  a  constant  amplitude.  The  electric  field  behavior  is 

similar,  but  the  peak  amplitude  in  this  case  is  much  smaller.  E  /(B  c)  = 

w  0 

-4 

1.4x10  .  In  a  one  dimensional  system,  the  peak  amplitude  is 

representative  of  saturation,  since  the  wave  energy  cannot  be  lost  by 
coupling  to  obliquely  propagating  waves  and  therefore  remains  constant. 
Throughout  the  simulation  run  less  than  0.2%  of  the  total  energy  is  lost 
demonstrating  that  computing  errors  are  negligible  and  that  the  energy  of 
the  system  is  conserved. 

The  time  history  of  the  different  particle  species  (Figure  6)  is  a 
most  enlightening  and  interesting  plot.  Panels  (a)  and  (b)  show  the 
behavior  of  the  hot  anisotropic  ions  as  a  functron  of  time.  In  panel  (a) 
we  see  that  the  hot  protons  are  losing  energy  in  a  direction  perpendicular 
to  the  ambient  magnetic  field.  This  results  in  wave  growth  which  in  turn 
causes  an  increase  in  the  hot  proton  parallel  energy.  This  expected 
tendency  towards  an  isotropic  hot  proton  distribution  is  best  seen  in  Panel 
(b)  where  we  plot  (T  as  a  function  of  time.  It  is  interesting  to 

note  that  the  temperature  anisotropy  decreases  rapidly  until  Q^t  =  360. 
reaching  an  anisotropy  of  (T =  1.3.  At  later  times,  the  anisotropy 
does  not  change  significantly,  but  reaches  a  value  of  1.2  at  the  end  of  the 
run.  Panels  (c)  and  (d)  show  the  large  heating  suffered  by  the  helium 
ions.  In  fact,  we  see  that  the  thermal  helium  has  been  heated  about  100 
times  Its  initial  value.  We  also  note  that  most  of  the  heating  :s  in  the 
perpendicular  direction.  Thermal  protons  are  also  heated  but  to  a  much 
lesser  degree. 

The  snapshots  of  phase  space  plots  v  -x.  v  -x  at  different  times  for 

X  y 

helium  ions  (Figure  7)  and  thermal  protons  (Figure  8)  show  complex 
nonlinear  structure.  Initially,  at  time  Q^t  =  0.  the  helium  ions  are 
uniformly  distributed  in  phase  space.  At  later  times.  Q.,t  = 


100.  we  see  a 


perturbation  in  v  .  ThiS  perturbation  occurs  from  a  simple  t  ,  x  force, 
and  results  in  a  transverse  motion  of  the  plasma.  Later  on.  Q^t  =  300  we 
see  that  the  v^  is  about  the  same  except  that  the  amplitude  of  the 
oscillations  is  larger,  and  tends  towards  smaller  wavenumbers.  At  that 
time,  we  note  osci nations  in  the  parallel  velocity  In  fact  the 
oscillations  appear  to  be  folding  on  themselves.  This  is  because,  as  the 
perturbation  in  the  transverse  velocity  occurs,  the  particles  experience  a 
force  6Vy  2  ^  y  parallel  direction.  At  C^t  =  400.  the  motion  is 
extremely  nonlinear,  with  large  striations  and  large  osci nations.  The 
motion  IS  so  nonlinear  that  a  mixing  of  plasma,  which  was  originally  in 
different  locations,  is  occurring,  resulting  in  a  hotter  distribution. 
Similar  nonlinear  motions  are  present  in  the  hydrogen  phase  space  plots 
(Figure  8)  but  they  are  not  as  prevalent  as  for  helium  ions  (Figure  7). 
The  reason  for  that  can  be  seen  from  the  linearized  momentum  equation  of 
motion  for  each  species. 


TT  '  *  “j  *  Bo' 


and  Maxwell  equation 


V  = 


Solving  for  v^,  we  get 


k  B  (1  -  u)/Q  ) 
o  j 

wf.ere  Q  =  (q/m)  B  . 

J  JO 

Thus,  the  species  response  is  proportional  to  B  and  depends  on  the 

w 

mass  through  the  resonant  denominator.  This  is  why  the  effect  is  less 

dominant  for  the  case  of  a  hydrogen  plasma.  As  the  plasma  evolves,  v  x  B 

J  w 

forces  become  substantial,  resulting  in  a  nonlinear  motion  along  the  field 


lines.  This  accounts  for  the  late  time  behavior  when  'spirals'  develop  in 
the  phase  space  and  to  a  lesser  extent  in  hydrogen. 

This  'weaker'  spiraling  for  thermal  protons  is  caused  by  the  fact  that 
the  fastest  growing  wave  is  in  the  LF  branch,  so  that  the  resonant  velocity 
is  far  from  the  thermal  velocity  of  protons.  In  Figure  9  we  show  the 
Fourier  analysis  of  the  wave  properties,  during  the  nonlinear  stage  of  the 
simulation  Qj^t  =  800  -1200.  We  have  decomposed  the  wave  into  the  forward 
and  backward  traveling  components.  It  is  found  that  both  the  forward  and 
backward  traveling  waves  have  equal  amplitudes,  and  only  the  LF  branch  has 
a  considerable  amplitude.  The  maximum  growth  rate  r/Q^  =  0.016  has  shifted 
to  lower  wavenumbers.  The  dominant  wavenumbers  are  best  seen  in  Figure  10. 
where  in  the  upper  panel  we  plot  a  three-dimensional  plot  of  the  amplitude 
of  the  wave  magnetic  field  versus  k.  as  a  function  of  time.  Looking  at 
Panel  (a)  it  is  easy  to  see  that  as  time  proceeds  the  wave  with  the  largest 
amplitude  occurs  at  longer  wavelength.  This  tendency  towards  longer 
wavelengths  can  also  be  seen  in  the  lower  panel  where  we  plot  the  time 
history  of  the  amplitudes  of  the  four  dominant  modes.  Initially,  mode  3 
and  4  seem  to  be  mostly  responsible  for  the  large  amplitude  wave  observed; 
however,  at  later  times  the  amplitude  of  mode  2  is  larger.  From  linear 
theory  analysis  the  shifting  of  the  fastest  growing  modes  to  lower 
wavenumbers  could  be  due  to  either  the  heating  of  the  helium  ions  or  to  the 
decrease  in  t  mperature  anisotropy  of  hot  protons  IGendrin  et  al . .  !984. 
Figure  (7a).  Figure  (7b)]. 

However,  as  can  be  seen  on  the  referenced  figures,  the  reduction  in 
the  bandwidth  of  the  unstable  waves  is  more  drastic  m  the  first  case  than 
in  the  second.  Besides,  since  the  most  amplified  wave  belong  to  the  LF 
branch  of  the  dispersion  curve,  a  decrease  of  the  hot  proton  anisotropy  has 
almost  no  consequence  on  the  frequency  of  the  maximum  growth  as  long  as  the 
anisotropy  remains  above  the  critical  value  [Kennel  and  Petschek.  1966] 
which  corresponds  to  it.  The  critical  anisotropy  is  defined  by: 


+  !  ' 


(11) 


where  =  T  /T  -1 . 
c  i 

For  (i)/Q^  =  O.I.  ~  0.2.  Consequently  the  frequency  of  the  fastest 

growing  wave  is  not  much  affected  by  a  decrease  in  ti.e  hot  proton 

anisotropy  as  long  as  (T  /T„),  >  1.2.  On  the  other  hand,  an  increase  of 

i  II  n 

the  hot  proton  parallel  energy,  as  evidenced  by  the  simulation  (see 
Figure  6).  may  be  at  the  origin  of  the  decrease  of  the  most  amplified  wave 
frequency  (see  e.g.  [Gendrin  et  al..  1971.  Figure  61). 

To  test  which  of  these  factors  is  responsible,  we  ran  the  simulation 
by  starting  with  hotter  thermal  proton  and  helium  distributions.  T,  =  T,,  = 

17  eV.  while  keeping  all  the  other  parameters  consent.  In  this  later 
simulation,  the  helium  ions  were  only  heated  by  a  factor  of  2.  while  the 
hot  proton  anisotropy  decreased  from  2  to  1.25  (still  larger  than  1.2)  and 
while  their  parallel  energy  Increased  by  a  factor  of  1.4.  The  shift  to 
lower  wavenumber  was  still  observed.  These  results  suggest  that  neither 
the  heating  of  the  helium  ions  nor  the  decrease  in  temperature  anisotropy 
of  the  hot  protons  is  responsible  for  the  lowering  of  the  most  amplified 
wavenumber.  We  therefore  conclude  that  the  shift  to  lower  wavenumber  is 
due  to  the  Increase  of  the  hot  proton  parallel  energy. 
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IV  Discussion  and  Conclusion 

In  this  paper  using  simulation,  we  have  studied  the  bulk  heating  of 
helium  ions  starting  from  an  initial  anisotropic  proton  distribution.  We 
have  shown  that  the  thermal  ions  are  heated  to  one  hundred  times  their 
initial  temperature,  due  to  the  growth  of  large  amplitude  ion  cyclotron 
waves.  Before  discussing  the  heating  mechanism  of  the  helium  ions  and  the 
saturation  of  the  ion  cyclotron  waves,  let  us  briefly  summarize  our 

findings. 

Wave  properties: 

1)  Large  amplitude  ion  cyclotron  waves  with  peak  amplitude  B  /B  =  0.06  are 

W  0 

observed  to  grow,  due  to  the  temperature  anisotropy  of  the  hot  proton 
distribution. 

2)  Analysis  of  dispersion  properties,  frequency  and  wavenumbers,  shows  that 

the  presence  of  helium  breaks  the  dispersion  relation  into  two 
branches,  one  below  the  helium  gyrof requency  (LF  branch)  and  one  above 
the  helium  gyrofrequency  (HF  branch).  This  is  in  agreement  with 
previous  work  on  the  linear  theory  dispersion  isee.  for  example. 

Gendrin  and  Roux  C19801). 

3)  At  the  early  stage  of  the  simulation,  agreement  with  linear  theory  is 

very  good.  The  fastest  growing  mode  occurs  on  the  branch  with  w-Uj 
^  0.16.  kV, 

=  0.17.  kV^/C^  =  0.31  and  r/C^  =0.022. 

4)  The  magnetic  field  reaches  a  peak  amplitude  at  rime  il.t  =  .50  and 

n 

oscillates  about  a  constant  value  thereafter. 

5)  At  later  times,  the  fastest  growing  mode  shifts  towards  longer 

wavelengths  but  remains  in  the  low  freque.ncy  branch.  'petra. 

t  =  300  and  t  =  1200  '.nows 

that  the  fastest  growing  mode  is  at  id/Q,,  =  1.1  and  kV,  Q,,  =  ts 

ri  AH 

discussed  in  section  111.  we  believe  this  shift  to  a  lower  n  mcae  .s 
due  to  an  increase  of  the  hot  proton  parallel  energy. 


analysis  for  the  time  period  between 


"h 


=  0.3  and 


r/c^  = 


=  0.02.  where  as  linear  theory  precicts 


Particle  properties: 


1)  The  hot  anisotropic  distribution  tends  towards  isotropy  with  time.  The 

temperature  anisotropy  changes  from  T  /T ^  =  2  at  Q^t  =0  to  T  /T ^  = 
1.3  at  Qj^t  =  360  and  then  changes  only  slightly  until  the  end  of  the 
run. 

2)  The  thermal  protons  are  only  slightly  heated  during  the  run.  As  noted 

previously,  this  is  because  of  the  fact  that  the  fastest  growing  wave 
IS  in  LF  branch,  causing  the  resonance  velocity  to  be  far  from  the 
thermal  velocity  of  the  protons. 

3)  The  cold  helium  ions  are  greatly  heated,  mostly  in  the  perpendicular 

direction.  The  perpendicular  helium  temperature  is  about  150  eV 
whereas  the  parallel  temperature  is  about  70  eV.  vtost  of  the 
perpendicular  heating  occurs  before  Q^t  =  360.  whereas  the  parallel 
heating  continues  until  Q^t  ^  600. 

4)  Phase  space  plots  for  helium  ions,  and  to  a  lesser  extent  hydrogen  ions. 

show  an  interesting  structure.  .At  the  early  stages  we  see  oscillation 
in  the  Vy  component,  which  is  simply  due  to  an  x  force.  As  time 
evolves  these  oscillations  grow  causing  a  significant  Sv^  ^  ^  y 
force  in  the  v.  direction.  This  results  in  spiraling  and  complex 
nonlinear  motion  in  the  v,|-x  phase  space  plots.  At  later  times,  we 
observe  thermal  nation  of  the  helium  ion  distributions. 

With  these  observations  in  mind  we  can  now  discuss  the  saturation 
mechanism  and  the  heating  of  helium  ions.  First  let  us  discuss  saturation 
of  the  wave.  If  we  calculate  the  linear  growth  rate  using  plasma 
parameters  at  the  time  of  peak  amplitude  Q^t  =  250.  we  find  a  positive 
growth  rate.  In  fact,  for  hot  proton  temperature  T„  =  20  keV.  T  /T„  =  1.5. 
cold  proton  temperature  T.  =  6  eV.  T  /T ^  =  2.  and  helium  temperature  T  ,  = 
30  eV.  T  'T  =  3.3.  we  find  from  linear  theory  calculation  r/Q^  =  0.012. 

The  growth  due  to  positive  anisotropy  must  be  balanced  by  a  damping 
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Tiecnanisrr.  One  obvious  candidate  for  this  i-Tping  ;s  trapping  of  heliums 
by  electromagnetic  waves.  Evidence  for  particle  trapping  can  be  seen  in 
the  phase  space  plot  v^-x  in  Figure  7.  In  fact  we  calculate  the  trapping 
frequency  for  the  mode  4  wave 


substituting  for  k  =  0.4,  B  =0.02,  q/m  =  0.25  and  v  =0.1.  at  Q„t  =  250. 

w  *  n 

we  find  =014.  Equating  ~  ^^arap  ’  damping  rate 

IS  of  the  same  order  of  magnitude  as  the  growth  rate. 

Now  we  turn  to  the  problem  of  heating  the  helium  ions.  First,  before 
discussing  this,  we  would  like  to  find  out  whether  the  heating  is  bulk 
heating  or  simply  high  energy  tail  formation.  To  do  that  we  plot  in 
Figure  II  the  distribution  functions  of  helium  ions  in  the  parallel 
direction  as  well  as  the  distribution  of  helium  ions  :r.  v^,,  one  of  the 
perpendicular  components,  at  various  times.  The  dashed  lines  are  the 
i.nitial  distribution  function.  Looking  at  the  left  column,  we  .note  that  at 
time  Q^t  =  20  there  is  some  heating.  However,  at  later  times  the 

distribution  is  much  hotter,  vioreover.  this  plot  shows  that  inaeed  we  do 
have  bulk  heating,  as  opposed  to  high  energy  tail  formation. 

Now  looking  at  the  time  history  plot  (Figure  6)  we  see  that  the 
heating  of  helium  seems  to  start  from  the  very  beginning.  Yet  the  thermal 
velocity  of  helium  is  such  that  it  is  difficult  to  understand  how  heating 
takes  place.  To  understand  the  exchange  of  energy  between  the  helium  ions 
and  the  wave,  we  calculate  E  "V^.  In  the  regions  where  E  -v  is  positive, 
the  helium  ions  are  gaming  energy  from  the  waves  and  v.ci  versa.  We  note 
that  in  general  during  the  initial  phase.  E  -v  is  positive  whereas  it  is 
generally  negative  later.  What  we  believe  is  happening  is  that  during  the 
initial  phase,  while  particles  are  accelerated  by  the  wave  there  is  no 
heating  since  the  nelium  ions  are  far  from  the  resonance  velocity.  We 


herefore  advance  the  following  scenario.  Initially,  the  helium  ions  are 


set  into  oscillation  (no  temperature  increase)  in  noth  the  perpendicular 
and  parallel  directions  by  fields  of  the  growing  wave.  When  v  ,  reaches  a 
value  near  the  resonance  velocity  such  that  heating  can  occur,  we  should 
see  a  sudden  increase  in  the  perpendicular  temperature.  We  define  a 
resonance  region  in  the  velocity  phase  space  as 


V'o  -  <  V  <  V’  +  V 

R  trap  i  R  trap 


(13) 


where  V’  is  the  resonance  velocity  given  by  (u  -Q,,  )/k.  and  V.  ^  is  the 
R  He  trap 

trapping  velocity  given  by  _ /k.  In  the  resonance  region  nonlinear 

trap 

trapping  of  heliums  is  possible,  leading  to  the  heating.  The  above 
scenario  is  consistent  with  the  lower  panel  of  Figure  6.  where  the  cold 
helium  temperature  increases  from  the  beginning  of  the  run.  The 
temperature  plotted  there  is  actually  the  space-averaged  temperature,  and 
the  oscillatory  effects  of  the  waves  have  been  averaged  over  many 
wavelengths  giving  an  apparent  velocity  spread  and  temperature  in  a  region 
of  space  a  quarter  of  the  wavelength  of  mode  4  in  each  direction,  in 
Figure  12  we  show  such  a  local  perpendicular  temperature  diagnostic  and 
also  the  number  of  helium  ions  whose  velocity  lies  within  the  resonance 
region  for  the  mode  4  wave  =  -0.16,  ^  0.07).  At  first  there  are 

no  particles  and  there  is  no  increase  in  temperature.  Later  the  number  of 
particles  has  increased  and  we  have  a  sudden  surge  in  temperature.  These 
diagnostics  clearly  indicate  that  the  energy  gained  by  the  heliums  during 
the  linear  phase  is  transfered  to  thermal  motion  due  to  the  nonlinear 
trapping. 

The  results  which  have  been  presented  show  the  efficiency  of  the 
simulation  code  to  interpret  complex  wave  particle  interactions  in  a 
multicomponent  plasma.  Some  interesting  results  have  been  obtained  which 
concern  the  heating  of  heavy  ions  in  the  equatorial  magnetosphere  by 
IGi-''s.  However,  more  e.xpen mental  runs  have  to  be  performed  ana  some 


improvements  are  still  needed  to  answer  specific  questions  related  either 


with  the  physics  of  the  phenomenon  or  to  the  efficiency  of  the  simulation 
code. 

As  far  as  physics  is  concerned,  there  remain  some  discrepancies 
between  the  simulation  results  and  experimental  data.  A  drastic  lowering 
in  the  frequency  of  the  most  amplified  wave  found  in  the  simulation  is 
generally  not  observed  in  experimental  data.  The  delay  which  is  observed 
between  the  establishment  of  a  rather  large  wave  field  and  the  heating  of 
He*  ion  IS  much  shorter  in  the  simulation  run  than  it  seems  to  be  in  the 
case  of  experimental  data  [Young  et  al . .  1981:  Roux  et  al..  19821.  The 
bunching  of  He*  ions,  which  is  a  striking  feature  of  ATS-6  observations 
[Mauk  et  al . .  19811  is  not  as  clearly  reproduced  in  the  results  of  the 
simulation. 

Some  of  these  discrepancies  can  be  explained  by  the  fact  that  the 
numerical  runs  reproduce  a  stationary  case,  and  not  a  progressive  one. 
because  a  periodic  boundary  condition  is  assumed  in  the  present 
simulation.  Waves  cannot  escape,  and  new  hot  particles  are  not  injected, 
thus  preventing  the  simulation  system  to  represent  an  equilibrium  where 
input  and  output  energy  flows  are  balanced.  The  size  of  the  system  should 
be  increased  (L^  >  128)  in  order  to  be  able  to  follow  in  more  detail  the 
evolution  of  the  most  amplified  wavenumber.  Finally,  a  parametric  study 
(number  of  particles  per  cell,  more  efficient  predictor-corrector  schemes. 
...)  should  be  undertaken  in  order  to  optimize  the  cost/eff iciency  of  the 
operational  code  which  has  been  so  successfully  used  throughout  this 


study. 


Figure  1.  Time  historv  of  the  magnetic  and  electric  energy  densities 
during  the  early  stages  of  the  simulation.  We  note  that  both  the 
magnetic  field  and  electric  fields  are  still  growing  at  time  IT^t  = 
200.  The  dashed  line  on  the  lower  panel  reflects  the  variation  of  E^.. 
the  smallness  of  which  demonstrates  the  efficiency  of  the  code  and  the 
electromagnetic  nature  of  the  wave. 

Figure  2.  The  time  history  of  the  energy  densities  of  the  hot  protons 

(panel  a),  helium  ions  (panel  b)  and  thermal  hydrogen  i panel  c).  ''iote 

that  the  temperature  anisotropy  hardly  changes  until  time  2.,i  =  IQO 

.d 

and  then  the  decrease  is  more  rapid.  Both  thermal  species  show  an 
increase  in  energy  density,  preferentially  in  the  perpendicular 
direction. 

Figure  3.  We  plot  the  dispersion  relation  of  the  component  during  the 
early  phase  of  the  simulation  a,_^t  =  204.  We  also  plot  the  theoretical 
linear  dispersion  relation  by  dashed  lines  using  the  Plasma  parameters 
at  Qj^t  =  100.  We  observe  that  both  the  simulation  results  and  the 
linear  theory  show  the  dispersion  relation  separated  into  two 
branches:  the  high  frequency  branch  (HF)  above  the  helium  cyclotron 
frequency  and  the  low  frequency  branch  (LF)  below  the  helium  cyclotron 
frequ  -".cy.  The  fastest  growing  frequency  occurs  in  the  LF  branch  at 
(u/Q^  =  0.16  and  =0.29. 

Figure  4.  In  the  upper  panel  we  plot  time  evolution  of  the  four  dominant 
modes  during  the  linear  phase  of  the  simulation.  The  lower  panel 
shows  the  theoretical  linear  growth  rate  as  a  function  of  wavenumber. 
From  the  upper  panel  we  find  that  the  growth  rates  for  mode  3  and  4. 
corresponding  to  'Q^  =  0.3  and  kV,/Ci^  =0.4  respect  i  vely,  are  r  = 
0.02.  The  lower  panel  predicts  the  peak  growth  rate  tc  be  r 
0.022  at  kV^/Q.^  =0,3. 

Figure  5.  Time  history  of  the  electric  and  magnetic  energy  densities  unt.l 


the  end  of  the  run.  The  magnetic  energy  aensity  increases  unt;!  about 

Out  =  300.  reaching  a  peak  amplitude  B  B  =0.06.  Saturation  occurs 

H  wo 

about  G^t  =  300, 

Figure  6.  .A  time  history  of  the  different  particle  species,  the  hot 

protons  panels  la)  and  (b).  the  perpendicular  temperature  of  helium 
and  thermal  protons  panel  (c)  and  the  parallel  temperature  of  helium 
and  thermal  protons  panel  (d).  We  note  in  panels  (a)  and  ib)  that  the 
anisotropy  decreases  rapidly  until  about  Q^t  =  360.  reaching  a  value 
of  T  =  1.3.  At  later  times  the  He^  ions  have  been  heated  about 
100  times  the  initial  value  in  the  perpendicular  direction.  Wost  of 
the  heating  occurs  prior  to  Qj^t  =  500.  The  ions  are  also  heated  in 
the  parallel  but  to  a  lesser  degree.  Thermal  protons  are  only 

slightly  heated. 

Figure  7.  Phase  space  plot  for  He'*’  at  different  times.  The  left  column  is 
:n  v'.-.x  space,  while  the  right  column  is  in  v  -x  space. 

>  X 

Figure  8.  A  phase  space  plot  for  thermal  hydrogen  at  different  times.  The 
left  column  is  in  v^-x  space  while  the  right  column  is  in  v^-x  space. 

Figure  9.  A  figure  of  the  dispersion  relation  of  the  B^  component  during 
the  nonlinear  stage  of  the  simulation  Q^t  =  800-1200.  The  upper  panel 
IS  for  the  forward  traveling  wave,  while  the  lower  panel  :s  for  the 
backward  traveling  wave.  At  this  stage  of  the  simulation  only  the  LF 
branch  has  a  considerable  stimulation.  The  dashed  line  represents  the 
linear  dispersion  relation  calculated  from  the  parameters  at  Q^t 
=1000.  Figure  iO.  Panel  (a)  shows  the  time  evolution  of  the  k- 
spectra  of  the  wave  magnetic  field,  while  panel  fb)  shows  the 
amplitude  versus  time  of  the  dominant  modes  in  the  system.  Sote  how 
the  fastest  growing  modes  shift  to  lower  wavenumbers  with  time. 

“igure  11.  A  plot  of  the  d!stri.oution  function  of  thermal  helium  at  three 
different  times.  The  left  column  is  the  parallel  distribution 
function,  whereas  the  right  column  is  the  distribution  function  in  one 


■  y.y. 
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of  the  transverse  directions.  The  dashed  line  is  the  initial 
distribution  function  at  time  Q^^t  =  0.  Note  the  broadening  of  the 
distribution  function  fiv^)  at  time  Qj^t  =  400.  Looking  at  the  lower 
left  panel  we  see  that  the  helium  continues  to  heat  until  the  end  of 
the  run.  but  at  a  slower  rate. 

Figure  12.  Here  we  plot  certain  diagnostics  needed  to  understand  the 
heating  mechanism  of  thermal  helium.  At  panel  (a),  we  plot  E  -v  . 
For  E  -v  >0,  the  energy  is  transferred  from  the  waves  to  the  helium 

X  i 

ions  and  vice  versa.  We  note  that  most  of  the  energy  transfer  from 
the  waves  to  the  helium  ions  occurs  prior  to  time  Q^^t  =  400.  In  panel 
(bi  we  plot  the  effective  temperature  of  the  helium,  by  measuring  the 
temperature  in  a  region  of  space  a  quarter  of  he  wavelength  of  mode  4 
in  each  direction.  In  panel  (c)  we  plot  the  number  of  helium  ions 
whose  velocity  lies  within  the  resonance  region.  We  note  that  number 
of  trapped  particles  and  the  effective  temperature  seem  to  increase 
simultaneously,  at  time  around  Qj^t  =  500. 
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igure  ..  Time  history  of  :;he  magnetic  and  electric  energy  densities  during  the 
early  stages  of  the  simulation.  We  note  that  both  the  magnetic  field  and 
ele^..rj.c  fields  are  still  growing  at  time  .‘jjt  ■  20C .  The  dashed  line  on  the 
lower  panel  reflects  the  variation  of  smallness  of  which  demonstrates 

t.oe  erriciency  of  :he  code  and  the  electromagnetic  nature  of  the  wave. 
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Figure  2.  The  time  history  of  the  energy  densities  of  the  hot  protons  (panel 
a),  helium  ions  (panel  b)  and  thermal  hydrogen  (panel  c).  Note  that  the 
temperature  anisotropy  hardly  changes  until  time  =  100  and  then  the  decrease 

is  more  rapid.  3oth  thermal  species  show  an  increase  in  energy  density, 
nref erentially  in  the  perpendicular  direction. 
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Figure  3.  We  plot  the  dispersion  relation  of  the  3y  component  during  the  early 
phase  or  the  simulation  iii^t  *  204.  We  also  plot  the  theoretical  linear 
dispersion  relation  by  dashed  lines  using  the  plasma  parameters  at  Fi^t  •  100. 

We  observe  that  both  the  simulation  results  and  the  linear  theory  show  the 
dispersion  relation  separated  into  two  branches;  the  high  frequency  branch  (HF) 
above  the  helium  cyclotron  frequency  and  the  low  frequency  branch  (LF)  below  the 
helium  cyclotron  frequency.  The  fastest  growing  frequency  occurs  in  the  LF 
branch  at  »  0.16  and  kVi/Slu  »0.29. 


(a)  Time  History  of  the  Wave  Magnetic  Field 
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(b)  Linear  Growth  Rates 
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Figure  4.  (a)  Time  evolution  of  the  four  dominant  modes  during  the  Linear  phase 

of  the  simulation,  (b)  Theoretical  linear  growth  rate  as  a  function  of  wave 
number.  From  Figure  4a  we  find  that  the  growth  rates  for  modes  3  and  4,  corres¬ 
ponding  to  “0.3  and  =  0.4,  respectively,  are  both  equal  to  = 

0.02.  Figure  4b  predicts  the  peak  growth  rate  to  be  Y/‘iH  “  0.022  at  kV^/ilH  “ 
0.3. 
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Magnetic  Energy  Density 
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Figure  5.  Time  history  of  the  electric  and  magnetic  energy  densities  until  th 
end  of  Che  run.  The  magnetic  energy  density  increases  until  about  H^t  *  300, 
reaching  a  peak  amplitude  By/Bo  -0.06.  Saturation  occurs  about  ilj^t  »  300. 
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Figure  6.  A  time  history  of  the  different  particle  species,  the  hoc  protons 
panels  (a)  and  (b),  the  perpendicular  temperature  of  helium  and  thermal  protons 
panel  (c)  and  the  parallel  temperature  of  helium  and  thermal  protons  panel  (d). 
Ve  note  in  panels  (a)  and  (b)  that  the  anisotropy  decreases  rapidly  until  about 
■l^t  “  360,  reaching  a  value  of  T_/T:i  •  1.3.  At  later  times  the  He*  ions  have 
been  heated  about  100  times  the  initial  value  in  the  perpendicular  direction. 
Most  of  the  heating  occurs  prior  to  ."li^t  =  500.  The  ions  are  also  heated  in 
the  parallel  but  to  a  lesser  degree.  Thermal  protons  are  only  slightly  heated. 
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Figure  9.  figure  or  the  dispersion  relation  of  the  By  conponent  during  the 
nonlinear  stage  of  the  simulation  rti|t  =  300-L200.  The  upper  panel  is  for  the 
forward  traveling  wave,  while  the  lower  panel  is  for  the  backward  traveling 
wave.  .\t  this  stage  of  the  simulation  only  the  LF  branch  has  a  considerable 
stimulation.  The  dashed  line  represents  the  linear  dispersion  relation 
calculated  from  the  parameters  at  ilut  -1000. 
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Figure  11.  A  plot  of  the  distribution  function  of  thermal  helium  at  three 
different  times.  The  left  column  is  the  parallel  distribution  function,  whereas 
the  right  column  is  the  distribution  function  in  one  of  the  transverse 
directions.  The  dashed  line  is  the  initial  distribution  function  at  time  ,;-j^t  = 
0.  Note  the  broadening  of  the  distribution  function  f(v  )  at  time  o^t  •  400. 
Looking  at  the  lower  left  panel  we  see  that  the  helium  continues  to  heat  until 
the  end  of  the  run,  but  at  a  slower  rate. 
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Figure  12.  Here  we  plot  certain  diagnostics  needed  to  understand  tne  heating 
mechanism  of  thermal  helium.  At  panel  (a),  we  plot  For  >  0,  the 

energy  is  transferred  from  the  waves  to  the  helium  ions  and  vice  versa,  'v'e  note 
that  most  of  the  energy  transfer  from  the  waves  to  the  helium  ions  occurs  prior 
to  time  riut  =  400.  In  panel  (b)  we  plot  the  effective  temperature  of  the 
helium,  by  measuring  the  temperature  in  a  region  of  space  a  quarter  of  ne 
wavelength  of  mode  4  in  each  direction.  In  panel  (c)  we  plot  the  number  of 
helium  ions  whose  velocity  lies  within  the  resonance  region.  We  note  that 
number  of  trapped  particles  and  the  effective  temperafire  seem  to  increase 
simultaneously,  at  time  around  ..ut  =  500. 
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